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Abstract 


We  solve  two  hydrodynamical  problems o  The  first  involves  a  shock 
wave,  a  contact  discontinuity,  and  a  rarefaction  wave  using  an  uncondi¬ 
tionally  stable  finite  difference  scheme.  The  Courant  condition  is  sat¬ 
isfied  everywhere  except  in  one  zone  behind  the  shock,  where  it  is  vio¬ 
lated  by  factors  of  10  and  100.  The  nonlinear  difference  equations  are 
solved  by  Newton's  method.  The  total  number  of  Newton  iterations  to  get 
to  a  certain  time  is  apparently  independent  of  the  degree  to  which  the 
normal  stability  condition  is  violated  in  the  one  zone. 

The  second  problem  involves  two  rarefaction  waves  moving  in  oppo¬ 
site  directions.  One  wave  moves  in  a  region  where  the  Courant  condition 
is  violated  by  a  factor  of  approximately  two.  The  other  wave  moves  in 
a  region  where  the  Courant  condition  is  satisfied.  Numerical  results 
are  compared  with  the  analytical  solution. 

An  examination  of  several  runs  indicates  one  explicit  time  step  is 
about  five  times  as  fast  as  one  implicit  time  step.  Therefore,  use  of 
the  implicit  method  is  indicated  when  the  Courant  condition  is  violated 


by  a  factor  of  5  or  more 
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Chapter  I 


Introduction 

Consider  a  hydrodynamieal  problem  in  which  a  shock  wave  or  dis¬ 
turbance  of  some  kind  is  advancing  into  a  material,.  Suppose  that  in 
the  neighborhood  of  the  disturbance  the  sound  speed  is  C^  and  suppose 
also  that  there  is  a  relatively  quiescent  region  behind  the  disturbance 
in  which  the  sound  speed  is  C^.  Any  explicit  finite  difference  method 
will  require 

max  C^<  1 
Ax 

so  that  if  C^  »  Cq  one  will  be  forced  to  follow  the  uninteresting  de¬ 
tails  of  the  motion  in  the  quiescent  region.  An  unconditionally  stable 
finite  difference  method  would  be  useful  in  such  a  problem.  We  present 
such  a  method  for  the  equations  of  nonviscous  compressible  flow  in  one¬ 
dimensional  Lagrangian  coordinates. 

The  Differential  Equations 

The  Lagrangian  hydrodynamic  equations  with  time  t  and  mass  m  as 


independent  variables  are: 


"  3m  =  ®  (niass  equation) 


bn 

St 


+ 


(momentum  equation) 


Si 

St 


du 

p3£ 


=  o 


(energy  equation) 


The  dependent  variables  are: 

p  =  density 

u  =  velocity 

p  =  pressure 

I  =  internal  energy 

The  velocity  is  defined  by 


u 


dx 
=  St 


vhere  x  is  the  coordinate  of  an  element  of  fluid  in  the  laboratory 
frame.  Differentiating  this  velocity  equation  with  respect  to  mass 
we  see  from  the  mass  equation  that 

1  _  bx 
p  -  5m 


The  Difference  Equations 

To  form  difference  equations  from  the  differential  equations  we 

imagine  the  fluid  partitioned  into  cells  of  mass  m  where  j  =  1,  2,  . .  .,  J, 

J 

J  being  the  total  number  of  cells.  Subscripts  on  field  variables  denote 


the  value  of  that  particular  variable  at  that  space  point.  For  example, 


Vl/2 


denotes  the  right-hand  cell  boundary  velocity  of  the  j  cell. 


Superscripts  denote  time  steps.  For  example,  denotes  the 

temal  energy  of  the  j  cell  at  time  t  =  (n+l)At. 

We  make  the  following  difference  approximations: 


n+1  n 

Uj+l/2  “  Uj+] 


m.  +  m,  - 
3  j+1 


w)  p?  -  A 


m.  +  m.. 
3  J+l 


dx 

5t  =  u 


n+l 

x.i+l/ 


"  0UJ+l/2  +  (W5"j+V2 


1  _  dx 
o  5m 


n+i  _  n+1 

J+l/2  j-1/2 


5u 

■  p55 


I?*-1  -  f?  Spr 

1  1  *  * 


/  n+1  n+1  \ 

(Uj-l/2  ■  Vl/2j  + 


(1-0  )P? 


l>l/2  -  U j+l/ 2 


v 


11- 


» 


* 


This  form  of  the  difference  equations  was  chosen  because  it  gives 
a  fairly  simple  form  to  the  Jacobian  matrix.  We  expect  that  the  Newton 
iterative  method  would  work  just  as  well  for  other  forms  of  the  equa¬ 
tions. 

The  polytropic  gas  equation  of  state  is  used.  Also  a  pseudo¬ 
viscosity  term  is  added  to  the  pressure  to  spread  the  shock  front.  The 
pressure  term  then  takes  the  form  [l]: 


if 


n+1 


=  (r-Dp“+1i"+1 

J  J 


Here  y  is  a  constant  characteristic  of  the  gas  and  \  is  a  constant 
whose  choice  will  be  discussed  later. 


(6) 


Rewriting  Equations  (1)  and  (4)  we  have: 

4r  -  g(1-9)fe  - 


n+l 

Vl/2 


-  u 


j+l/  2 


m.  +  m.  _ 
0  0+1 


m.  +  m.  i 
j  J+l 


=  0 


i 
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V 


V 


(7)  I*?*1  -  I?  - 

0  J 


for  j  —  lj  2,  •  •  •  y  J. 

Assuming  that  ve  know  the  values  of  the  dependent  variables  at  time 
n,  this  gives  us  a  system  of  2J  simultaneous  nonlinear  equations  in  2J 
unknowns  for  the  values  at  time  n+1. 

Newton's  method  can  be  used  to  solve  this  system  of  equations.  For 
a  general  system  of  the  form 

f(y)  =  0 


where  f  and  y  are  vectors  Newton's  method  is  an  iterative  procedure  in 
which  the  p+l-st  iterative  is  defined  by 


y<^>  =  y<*> 


+  Ay 


where  Ay  is  the  solution  of  the  linear  system 
(8)  JAy=-f[y^p)] 
where 

J  = 

evaluated  at  y^p^. 

Taking  (5)  into  account  we  see  that  (6)  and  (j)  may  be  written  in 
the  form 


v 


•13' 


(9)  «j  (""W  ^  Ci)  -  ° 


In  oar  case  we  see  then  that  the  Jacobian  matrix  has  a  particularly 
simple  form,  namely  that  it  is  block  tridiagonal. 


where  the  submatrices  are  2x2. 


V 


v 


where  all  differentiations  are  with  respect  to  the  variables  u  or  I  at 
time  n+lo 

We  use  the  usual  scheme  to  invert  a  block  tridiagonal  matrix  [2]. 
Define  2x2  matrices  as  follows: 


W1  ■=  B-1^ 


-  A.W. 
3  J 


j-i)  cj’ 


2  <  j  <  J-! 


Gi  ■ 


(B.  -  A.W.  (f .  -  A.G. 

{  3  3  J“l]  [  3  3  J-l/J 


2  <  j  <  J 


If  we  redefine  A V  and  f  so  that 

fj  -  (gj>5o)>  ’  rJ+i/2'aij) 

then 


-15 


The  p+l~st  iterate  is  obtained  by  setting 


u(l*l)  =  U(P) 

Vl/2  Uj+l/2 


+  All 


j+1/2 


I*!*1* 

0 


+  AX . 
0 


Chapter  II 


Derivation  of  Matrix  Elements  and  Stability  Analysis 


We  now  derive  the  entries  of  the  2x2  submatrices  of  the  Jacobian 
matrix.  Referring  to  (6),  (7),  (ll),  (12),  and  (13)  we  have: 


V 
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V 


V 


-20  At 

mJ  +  “HI 


n  n+1 


*n+l 
*PH1 


H3/2  ^+3/2; 


-39  At 

"j  +  Vl 


'dp"1 

-LJ—  -  ■  V1 


0 


0 


Here  again  all  differentiations  with  respect  to  u  and  I  are  to  be  taken 
at  time  n+1* 

To  complete  the  derivation  we  need  the  various  partial  derivatives 
of  the  pressure  terms. 

From  (2)  it  follows  that 


,,  \  n+1  n+1  n 

l11*)  /O  -  X,_T  /O  =  *4, 


n 


‘,3+1/2  "  Xj-l/2  “  XHl/2  "  XHl/2  +  9At(Uj  ,n  "  u"+’ 


Hl/2  "  j-l/2j 


(l-0)At^+1/2  -  uj-ay2) 


The  pressure  term  may  then  be  written 


n+1  =  (y-1)^1^  /  79  3  (  n+1  n+1  \ 

PJ  xn+1.  .  xn+l  /n+1  _  n+1  (Uj-l/2  “  V/2 

Hl/2  xj-l/2  V  Hl/2  Xj-l/2 


if 


W2  -  UHV2  >  0 
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J 


and 


n+1 

P0 


— _ Ljla _ 

n+1  n+1 

Xj+l/2  "  Xj-l/2 


if 


n+1 

Uj-l/2 


n+1  ^ 

0+1/2  < 


A  tabulation  of  the  pressure  derivatives  follows: 


dp' 


n+1 


0+1/2 


(7-l)m  9AtI  /  Tra  p 

- s2 - H — _  -  x  / - SUL - 

(X,j+V2  -  *3- 1/2)  V  XJ+1/2  '  Xj-V2 


-X0At 


jSl 


0+1/2  “  Xo-l/2 


(“3-yg  -  Vi/a) 

(*3+1/2  "  *3-0/2) 


where  the  last  two  terms  do  not  appear  if  uj_i./2  “  uj+i/2  — 


dp 


n+1 

itL.  *  « 


dp' 


^+1/2 


n+1 

0+1 


^7+3/2 


dp^+1  (r-l)m.  Spj+1 

dI0  (X0+l/2  “  X0-l/2)  *  j 


^j+l/2  L~J 


n+1 
P  • 


dpn+1 

(Uj-l/2  "  Uj+l/2)  =  (Uj-l/2  "  Uj+l/2)  du^l/2  ’  Pj+1 


5 

3T 

0 


n+1/ n+1  n+1  \ 

PJ  (j-1/2  ’  Uj+l/2) 


dp 


n+1 


(Uj-l/2  "  uJ+l/2)  -it 


>  n+1  v  n+1  v  n 

dpo  Spi  dpi+l 

^0+1/2  ^0-1/2 


n+1. 


n  n+1 


3^  LPJ  (U>V2  -  U>1/2)J  -  ^  +  (Uj-l/2  -  Vl/2)  5Tj^ 


^r1  c  <:  <t 

^3/2  ’  *V3/2  ' 


^=0;  ^ - i^V 

;H1  5Id+l  Xj+3/2  -  Xj+l/2 


Thus  if  ve  vrite 


K  =  0At 


Am  = 


K 


m.  +  m.  , 
i  0+1 


^0  U j-l/ 2  "  Uj+l/2 
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we  have: 


The  method  described  in  Chapter  I  is  used  to  invert  this  matrix. 


Stability  Analysis 

As  has  been  pointed  out  [3]  a  rigorous  stability  analysis  for  the 
hydrodynamical  difference  equations  has  not  been  carried  out.  This 
analysis  proceeds  in  the  same  manner  as  that  done  by  Fromm  [l]. 


-21- 


We  assume  that  the  field  variables  vary  slightly  from  their  true 
values. 

Let: 


0+1/2  =  uo(x  +  *0+1/2) ’ 

5  «  1 

PJ =  M1  +  £j)’ 

€  «  1 

ZS  *  M1  +  bi)’ 

8  «  1 

For  simplicity  assume  that  the  cell  masses  are  equal;  m.  =  m  =  pn5xn<> 

We  substitute  these  values  in  the  difference  equations  and  get  the 
equations  of  first  variation,  dropping  all  higher  order  terms* 


m,1  Pq8*o 


XJ+V2  '  =  pn  (1+£.. 

a  1  j 


n\  °V 


Then  (14)  becomes 


<«>  ^+1  -  S  - ^  (•&,  -  <&8)  ♦  M 


05tu 


(l~0 )&tur 


We  define  the  Courant  number  p  = 
speed. 

Then  (15)  becomes 


C05t 

6x« 


where  CQ  is  the  local  sound 


(lfi,  cf-  -  ^  (S^2  .  |~^2)  *  ^2  (|“.V2  -  |^+V2) 


•22' 


The  first  variation  of  the  energy  equation  (7)  is 
/.„s  _  /cn+l  _n\  06t  n+1  /.n+1  .n+1  \ 

(17)  Jo  ft  -  Bj)  -  —  “oft-Vs  •  W) 


)5t  n  f„n  .n 

Pj  Uo(5j-l/2  **  ^ j+l/2 


For  p.  we  substitute 
0 


(7-ib0i0(i  -  5“  +  e“)  ♦  %Vo(f“-y2  •  s”+va) 


Upon  simplification  (17)  becomes 


(18)  sVx  -  &1?  . 

J  J 


(r-i)PnuA5t 


6  (^-1/2  "  ^+1/2)  +  t1-0)  (^j-i/2  ‘  ^+1/2) 


Finally  we  get  the  first  variation  of  the  momentum  equation  by  sub¬ 
stitution  into  (6)  and  again  dropping  higher  order  terms. 

<w>  -  %  {t^Vo  [(^r1  -  «S$ +  (6“+1  - <&)] 


+  Xp0°0u0  (lj-l/2  “  5J+l/2)  “  (^+1/2  "  5j+3/2) 


^  f  -  Vo  [ft  "  ^+l)  *  ft  "  6^)] 

Xp0C0U0  ft-l/2  “  ^+1/2)  “  ft+l/2  "  ^+3/2)] 
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At  this  point  ve  do  the  usual  thing.  We  assume  that  the  coeffi¬ 
cients  are  constant  and  that  the  solution  of  these  three  equations  of 
first  variation  can  be  vritten  in  a  Fourier  series.  If  so,  then  each 
term  of  the  series  is  a  solution  and  ve  look  at  a  typical  term  to  see 
vhat  conditions  must  be  satisfied  to  make  it  a  solution. 

We  assume  that 

&1/2  -  ^(J+1/2^ 


€ 


n 

d 


ee 


iki  n 
r2 


6°  =  be^r” 


and  consider  only  the  special  case  r^  =  =  r^  =  r. 

Substitution  of  these  values  into  (16),  (l8),  and  (19)  yields  after 
simplification 

(rO+1-0 )  2i  sin  k/2  jiuf 


(l-r)e  - 


^0 


t  =  0 


2iC^p  sin  k/2  (r@+l-@) 


V 


e  +  [l  -  4Xp  sin2  k/2  (r0+l-0)-r]| 
2iCQp  sin  k/2  (r@+l-0) 


V 


6  =  0 


2i(7-l)u0p  sin  k/2  (r0+l-0) 


'0 


|  +  (l-r)6  =  0 
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For  this  system  of  homogeneous  linear  equations  to  have  a  non¬ 
trivial  solution  it  is  necessary  that  its  coefficient  matrix  he  singular. 


2i(j.u0  sin  k/2  (rS+1-0) 

C0 

1  -  r 

0 

2iC sin  k/2  (r@+l-0) 

p 

2iC  \x  sin  k/2  (rG+1-0) 

0 

V 

1  -  4Xn  sin  k/2  (r0+l-9)-r  - 

2i(y-l)u-(i  sin  k/2  (rO+1-0 ) 

_ Q _ 1 - 

v 

0 

co 

1  -  r 

Expanding  this  determinant  we  get: 

(l-r)[l*r-4A,p  sin2  k/2  (r0+l-0)j  +  4p2  sin2  k/2  (r0+l-0)2  =  0 

For  full  generality  at  this  point  we  would  have  to  study  the  roots 
of  this  quadratic  equation  for  arbitrary  0.  This  is  somewhat  difficult. 

p  O 

The  two  cases  of  most  importance  are  0  =  1/2  and  0  =  1.  Let  sin  k/2  =  t  . 
For  0  =  1/2  the  equation  reduces  to 

r2(l  +  2A.pT2  +  p2T2)  +  rfep2T2  -  2)  +  (l  -  2A.pT2  +  p2x2)  =  0 


1  +  2A.pT  +  p  t' 
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V 


Case  (l): 


If  X2t2 


<  1,  then  r  is  complex 


uW  ±  4uV(x  -  A2) 

(l  +  2X|iT2  +  p2t2)2 


1  -  2lpT2  +  |i2T2 
1  +  2XpT2  +  (i2T2 


<  1 


Case  (2):  If  X2t2  >  1,  then  r  is  real. 

To  have  r  <  1  we  need 

2X(it2  +  |i2r2  >  -  p2T2  +  2|aT  */x2t2  -  1  or  Xt  +  pT  >  +  n/x2t2  -  1 

But  Xt  >  \/x2 t2  -  1,  so  indeed  r  <  1.  The  proof  for  the  case  0  =  1  is 
similar. 

Notice  here  that  r  <  1,  independent  of  p,  the  Courant  number.  This 
shows  that  we  have  verified  a  necessary  condition  for  this  method  to  be 
unconditionally  stable,  namely,  for  solutions  of  the  equations  of  first 
variation  having  the  form  we  have  prescribed. 
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Chapter  III 


Discussion  of  Numerical  Results 

The  first  problem  used  to  test  the  scheme  is  the  one  used  in  [2]. 
Here  we  have  two  constant  states  separated  by  a  discontinuity.  The  con¬ 
figuration  at  40  cycles  is  a  shock  moving  with  speed  1.24,  a  contact 
discontinuity  at  the  point  of  initial  discontinuity,  and  a  rarefaction 
wave. 

Hie  initial  conditions  for  this  problem  are: 
y  =  1.4 

At  =  0.337  (the  Courant  value) 

J  =  50  (25  cells  to  the  left  of  the  interface  and  25  cells  to  the 
right). 


Material  on  left 


Material  on  right 


Tables  1,  2 ,  and  3  give  the  velocities,  densities,  and  internal 
energies  for  several  different  calculations. 

The  Lax-Wendroff  figures  refer  to  the  values  obtained  using  the 
scheme  of  Reference  4. 

Exact  refers  to  the  analytical  values. 

Explicit  refers  to  values  obtained  using  one  of  the  explicit  schemes 
of  Reference  1. 

Imp^  refers  to  calculations  done  with  all  50  cells  having  mass  one. 

To  test  numerically  the  unconditional  stability  of  the  implicit  dif¬ 
ference  scheme  a  thin  cell  having  the  same  density,  but  only  a  tenth  the 
mass  and  width  of  the  other  cells,  was  put  into  cell  20.  This  means 
that  the  Courant  condition  was  violated  there  by  a  factor  of  approxi¬ 
mately  ten.  Inrpg  refers  to  calculations  done  with  this  thin  cell. 

Inrp^  is  similar  to  Imp^,  the  only  difference  being  that  this  time 
cell  20  was  given  mass  and  width  one-hundredth  that  of  the  other  cells. 
Thus  the  Courant  condition  was  violated  by  a  factor  of  approximately 
one  hundred.  As  can  be  seen  from  the  results  for  Impg  and  Imp^,  no 
instability  appeared  in  the  calculation.  When  the  explicit  method  was 
run  with  a  thin  cell,  large  fluctuations  appeared  and  eventually  two 
cell  boundaries  crossed  near  the  thin  cell. 

Since  Newton’s  method  involves  evaluating  the  elements  of  a  large 
matrix  and  then  inverting  it,  another  method  for  solving  the  system  of 
simultaneous  nonlinear  equations  was  considered,  namely,  the  method  of 
nonlinear  overrelaxation  as  described  in  [5].  If  one  has  a  system  of 
k  algebraic  equations 


Xp  ....  Xj  .  0,  i  =  l,2,...,k 

each  having  one  continuous  derivative,  then  the  generalization  of  ordi¬ 
nary  overrelaxation  suggested  by  Lieberstein  for  the  nonlinear  system 
is 


etc,,  where  f^  =  df^/bx^  Here  superscripts  on  variables  denote  the 
n  iterate  and  n+l-st  iterate  and  o>  is  the  relaxation  parameter. 

It  was  hoped  that  this  method  would  be  faster  than  Newton's  method 
for  solving  the  system  of  nonlinear  equations.  As  Lieberstein  points 
out,  the  rate  of  convergence  of  this  method  depends  rather  critically  on 
the  choice  of  a>.  For  our  choice  of  oo  =  1  the  overrelaxation  method  was 
actually  slower  than  Newton's  method,  but  a  more  careful  study  of  how  to 
choose  to  in  an  optional  manner  would  probably  make  the  overrelaxation 
method  faster  than  Newton's  method. 

Figure  1  gives  the  velocity  profile  for  Imp^  superimposed  on  the 
exact  solution.  Figure  2  gives  the  density  profile  for  Inrp^  superim¬ 
posed  on  the  exact  solution. 
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Figure  3  gives  the  velocity  profile  near  the  shock  front  for 
0=1/2  and  three  different  values  of  X,  In  general  a  large  value  of 
X  gives  a  smoother  profile  near  the  shock  front  but  spreads  the  shock 
over  several  cells.  A  smaller  X  gives  a  sharper  shock  front  but  has 
more  oscillation.  Some  intermediate  value  of  X  gives  the  best  compro¬ 
mise  between  these  two  effects.  We  have  found  that  for  0  =  1/2  a  smaller 
X  can  be  used  than  for  the  explicit  case.  His  is  clear  from  Figure  3« 

To  test  the  explicit  case  we  needed  to  take  X  =  0.7  and  even  then  some 
oscillations  appeared  near  the  shock  front,  but  for  the  implicit  case 
X  =  0.5  gave  a  fairly  sharp  shock  front  with  practically  no  oscillation. 

Several  trials  were  run  with  0=1.  The  most  notable  differences 
in  the  results  are  that  (l)  they  are  somewhat  less  accurate  than  for 
0=l/2  and  (2)  it  was  found  that  the  pseudo-viscosity  term  was  un¬ 
necessary  and  X  =  0  gave  the  sharpest  shock  front  with  little  oscil¬ 
lation. 

The  reduced  accuracy  may  be  understood  when  one  considers  that  for 
0=l/2  all  the  differences  are  centered  and  the  truncation  errors  are 
of  order  (At)  .  For  any  other  choice  of  0  some  second-order  truncation 
error  is  present.  One  should  then  expect  more  accurate  results  for 
0=l/2  than  for  any  other  choice  of  0. 

The  result  that  X  =  0  is  the  best  choice  means  that  the  implicit 
scheme  itself  contributes  an  effective  viscosity  term  when  0=1. 

Table  4  gives  the  velocity  profile  for  0=1,  X  =  0  and  X  =  0.5. 
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The  total  number  of  Newton  iterations  required  to  get  to  t  =  13»^8 
is  approximately  120  for  Imp^,  Imp^,  and  Imp^  and  for  9  =  0.5  and  0  =  1. 
This  number  is  thus  apparently  independent  of  the  degree  to  which  the 
Courant  condition  is  violated  in  zone  20. 

The  convergence  criterion  used  required  that  the  maximum  percentage 
change  in  any  value  of  p  or  I  be  less  than  1$  on  one  Newton  iteration. 
This  generally  required  three  Newton  iterations  for  each  time  cycle. 

When  this  criterion  was  relaxed  to  the  point  of  requiring  the  maximum 
percentage  change  to  be  less  than  10$  the  final  results  were  changed  at 
most  by  a  unit  or  two  in  the  fourth  significant  digit.  For  this  cri¬ 
terion  only  two  Newton  iterations  were  required  for  each  time  cycle. 

When  the  stricter  convergence  criterion  is  used  timing  experiments 
have  indicated  that  the  explicit  method  is  approximately  five  times 
faster  per  time  step  than  the  implicit  scheme.  In  this  case  use  of  the 
implicit  scheme  is  indicated  when  the  Courant  condition  is  to  be  vio¬ 
lated  by  a  factor  of  5  or  more. 

With  the  less  stringent  convergence  criterion  the  explicit  method 
is  only  about  three  times  as  fast  as  the  implicit  method.  Thus  if  this 
convergence  criterion  gives  sufficient  accuracy,  and  in  practice  it  has, 
use  of  implicit  scheme  is  indicated  when  the  Courant  condition  is  to  be 
violated  by  a  factor  of  3  or  more. 

In  practice  the  rate  of  convergence  has  been  approximately  quad¬ 
ratic,  the  maximum  percentage  change  being  roughly  squared  each  time. 

Also  in  practice  the  Jacobian  matrix  has  proved  to  be  diagonally  dominant. 
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Thi6  partially  accounts  for  the  accuracy  of  the  inversion  routine  and 
the  rapidity  of  convergence  of  Newton’s  method. 

The  second  problem  involves  one  gas,  half  of  which  is  initially  at 
rest,  the  other  half  initially  moving  with  constant  velocity.  The  con¬ 
figuration  at  40  cycles  is  two  rarefaction  waves  moving  in  opposite 
directions  at  a  sound  speed  which  is  C  =  0.316. 

The  gas  initially  at  rest  is  divided  into  90  cells  of  mass  0.1. 

The  gas  which  is  initially  moving  is  divided  into  10  cells  of  mass  1.0. 
Other  initial  conditions  are: 

7  =  1.4  Q  =  0.5 
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I.  =  0.1786 

J 

Ij  =  0.1786 

Figure  4  gives  the  plot  of  the  density  at  time  t  =  20  from  the 
numerical  results  and  also  the  analytical  values.  It  can  be  seen  that 
no  instabilities  have  appeared  in  the  finely  divided  material  even 
though  the  Courant  condition  is  violated  by  a  factor  of  approximately 
two. 
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One  difficulty  apparent  from  the  graph  is  that  the  implicit  scheme 
seems  to  lag  behind  the  true  solution  in  the  finely  divided  region. 
Evidently  the  scheme  does  not  allow  signals  to  be  propagated  at  sound 
speed  in  such  a  finely  divided  material.  This  may  be  the  fault  of  the 
form  of  the  difference  equations,  since  they  are  nonconservative. 
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Table  1  Velocity 


Cell 

Lax 

Wendroff 

Explicit 

Exact 

Imp^ 

Ing?2 

^3 

1 

0.702 

O.698 

0.698 

0.702 

— - 

— 

2 

0.709 

0.699 

0.698 

0.707 

O.709 

O.709 

3 

0.725 

0.703 

0.698 

0.720 

0.719 

O.721 

4 

0.754 

0.716 

0.698 

0.747 

0.743 

0.746 

5 

0.600 

0.752 

0.698 

0.794 

0.788 

0.793 

6 

0.866 

0.826 

0.822 

O.865 

O.857 

0.864 

7 

0.9^6 

0.938 

0.984 

0.961 

O.95O 

0.960 

8 

1.045 

1.075 

1.130 

1.075 

I.063 

1.074 

9 

1.150 

1.223 

1.342 

1.203 

1.190 

1.202 

10 

1.259 

1.372 

1.453 

1-335 

1.322 

1.334 

11 

1.366 

1.512 

I.528 

1.463 

1.452 

1.463 

12 

1.463 

1.616 

I.528 

1.568 

1.611 

1.568 

13 

1.541 

1.613 

I.528 

1.611 

I.588 

1.612 

14 

1.589 

I.588 

I.528 

1.586 

1.563 

1.587 

15 

1.596 

1.563 

I.528 

I.561 

1.546 

1.562 

16 

1.566 

1.548 

I.528 

1.545 

1.537 

1.546 

17 

1.525 

1.538 

1.528 

1.537 

1-533 

1-538 

18 

I.508 

1-533 

I.528 

1-533 

1.531 

1.534 

19 

1.518 

1-530 

1.528 

1.531 

1-531 

1.532 

20 

1.534 

1.528 

1.528 

1.529 

1-530 

1.529 

1.530 

1-530 
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21 

1.529 

1.526 

1.528 

1.528 

1.528 

1.529 

22 

1.526 

1.525 

1.528 

1.527 

1.527 

1.527 

23 

1.528 

1.527 

1.528 

1.526 

1.526 

1.527 

24 

I.528 

1.532 

1.528 

1.526 

1.526 

1.526 

25 

I.528 

1.520 

1.528 

1.526 

1.526 

1.526 

26 

I.528 

1.526 

1.528 

1.525 

1.526 

1.526 

27 

I.528 

1.527 

1.528 

1.525 

1.525 

1.525 

28 

1.528 

1.528 

1.528 

1.525 

1.525 

1.526 

29 

1.528 

1.527 

1.528 

1.525 

1.526 

1.527 

30 

I.528 

1.525 

1.528 

1.526 

1.526 

1.527 

31 

1.528 

1.529 

1.528 

1.526 

1.526 

1.526 

32 

I.527 

1.532 

1.528 

1.526 

1.526 

I.526 

33 

I.528 

1.522 

1.528 

1.525 

1.525 

1.526 

3^ 

1.527 

1.530 

1.528 

1.526 

1.525 

I.526 

35 

1.527 

1.512 

1.528 

1.525 

1.525 

1.526 

36 

I.528 

1.531 

1.528 

1.528 

1.527 

1.529 

37 

1.533 

1.524 

1.528 

1.524 

1.524 

1.525 

38 

1.526 

1.529 

1.528 

1.528 

1.529 

I.528 

39 

1.519 

1.472 

1.528 

1.508 

1.508 

1.509 

ko 

1.576 

1.518 

1.528 

1.522 

1.523 

1.527 

kl 

1.546 

1.298 

1.528 

1.432 

1.432 

1.423 

k2 

0.850 

0.725 

0 

0.558 

0.558 

0.547 

43 

0.108 

0.136 

0 

0.114 

0.114 

0.111 

44 

0.006 

0.018 

0 

0.019 

0.019 

0.019 

45 

0.000 

0.002 

0 

0.003 

0.003 

0.003 
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Table  2  Density 


Lax 


Cell 

Wendroff 

Explicit 

Exact 

Imp! 

Imp2 

Impj 

1 

0.445 

0.445 

0.445 

0.445 

— 

— 

2 

0.444 

O.445 

0.445 

0.445 

0.445 

0.445 

3 

0.442 

0.445 

O.445 

0.444 

0.444 

0.444 

1+ 

O.438 

0.444 

0.445 

0.441 

0.441 

0.441 

5 

0.432 

0.440 

0.445 

0.436 

0.437 

0.436 

6 

0.423 

0.431 

0.426 

0.428 

0.429 

0.429 

7 

0.413 

0..418 

0.407 

0.4l8 

0.419 

0.418 

8 

0.401 

0.402 

0.388 

O.4o4 

0.406 

0.405 

9 

0.388 

0.384 

0.370 

0.389 

0.391 

0.390 

10 

0.376 

0.367 

0.350 

0.374 

0.375 

0.374 

11 

0.363 

0.351 

0.345 

0.358 

0.360 

0.359 

12 

0.352 

0.338 

0.345 

0.345 

0.347 

0.346 

13 

0.3*6 

0.334 

0.345 

0.337 

0.337 

0.337 

14 

0.338 

0.338 

0.345 

0.336 

0.336 

0.336 

15 

0.337 

0.339 

0.345 

0.339 

0.339 

0.339 

16 

0.341 

0.341 

0.345 

0.342 

0.341 

0.342 

17 

0.345 

0.343 

0.345 

0.343 

0.343 

0.343 

18 

0.347 

0.344 

0.345 

0.344 

0.344 

0.344 

19 

0. 346 

0.345 

0.345 

0.344 

0.344 

0.344 

20 

0.344 

0.345 

0.345 

0.344 

0.345 

0.345 

0.344 

0.345 

21 

0.344 

0.345 

0645 

0.344 

0.345 

0.345 
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22 

0.345 

0.345 

0.345 

0.345 

0.345 

0.345 

23 

0.344 

0.345 

0.345 

0.345 

0.345 

0.345 

24 

0.344 

0.345 

0.345 

0.345 

0.345 

0.345 

25 

0.346 

0.345 

O.345 

0.345 

0.345 

0.345 

26 

1.212 

1.170 

1.304 

1.213 

1.218 

1.170 

27 

1.287 

1.225 

1.304 

1.242 

1.247 

1.225 

28 

1.295 

1.259 

1.304 

1.264 

1.269 

1.259 

29 

1.292 

1.284 

1.304 

1.280 

1.283 

1.284 

30 

1.294 

1.289 

1.304 

1.290 

1.292 

1.289 

31 

1.295 

I.289 

1.304 

1.296 

1.296 

1.299 

32 

1.297 

1.286 

1.304 

1.301 

1.301 

1.304 

33 

1.300 

1.296 

1.304 

1.304 

1.303 

1.306 

34 

1.302 

1.294 

1.304 

1.305 

1-305 

1.306 

35 

1.304 

1.304 

1.304 

1.306 

1.306 

1.308 

36 

1.306 

1.294 

1.304 

1.303 

1.303 

1.304 

37 

1.307 

1.292 

1.304 

1.305 

1.305 

1.306 

38 

1.306 

1.286 

1.304 

1.294 

1.294 

1.296 

39 

1.299 

1.307 

1.304 

1.306 

1.306 

1.309 

4o 

1-335 

1.248 

1.304 

1.261 

1.261 

1.266 

4l 

1.326 

1.194 

1.304 

1.313 

1-313 

1.312 

42 

0.831 

0.964 

0.500 

0.841 

0.841 

0.834 

43 

0.540 

0.628 

0.500 

O.568 

O.568 

O.566 

44 

0.503 

0.518 

0.500 

0.511 

O.512 

0.511 

45 

0.500 

0.502 

0.500 

0.502 

0.502 

0.502 
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Table  3  Internal  Energy 


Lax 


Cell 

Wendroff 

Explicit 

Imp-L 

Inrpg 

^3 

1 

19.78 

19.76 

19.76 

— 

— 

2 

19.77 

19.76 

19.75 

19.76 

19.76 

3 

19.74 

19.75 

19.72 

19.73 

19.73 

4 

19.67 

19.73 

19.68 

19.69 

19.68 

5 

19.55 

19.68 

19-59 

19.60 

19.59 

6 

19.40 

19.5^ 

19.45 

19.47 

19.45 

7 

19.21 

19.32 

19.26 

19.28 

19.26 

8 

18.98 

19*02 

19.01 

19.03 

19.01 

9 

18.74 

I8.69 

18.72 

18.75 

18.73 

10 

18.49 

18.34 

18.42 

18.45 

18.43 

11 

18.25 

18.01 

18.12 

18.15 

18.12 

12 

18.02 

17.73 

17.85 

17.87 

17.85 

13 

17.92 

17.60 

17.67 

17.68 

17.67 

14 

17.73 

17.65 

17.66 

17.65 

17.65 

15 

17.73 

17.70 

17.72 

17.71 

17.72 

16 

17.79 

17.75 

17.77 

17.77 

17.77 

17 

17.89 

17.79 

17.80 

17.80 

17.80 

18 

17.92 

17.81 

17.82 

17.82 

17.81 

19 

17.90 

17.82 

17.82 

17.82 

17.82 

20 

17.86 

17.82 

17.83 

17.83 

17.83 

17.83 

17.83 

21 

17.86 

17.82 

17.83 

17.83 

17.83 

22 

17.88 

17.83 

17.84 

17.83 

17.83 
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23 

17*90 

17.83 

17.84 

17.84 

17.84 

24 

17.90 

17.81 

17.84 

17-84 

17.84 

25 

17.78 

17.80 

17.86 

17.86 

17.85 

2 6 

5.086 

5.266 

5.074 

5.074 

5.052 

27 

4.789 

5.026 

4.954 

4.955 

4.938 

28 

4.763 

4.880 

4.871 

4.871 

4.852 

29 

4.770 

4.803 

4.813 

4.813 

4.797 

30 

4.765 

4.772 

4.772 

4.772 

4.761 

31 

4.757 

4.747 

4.747 

4.747 

4.738 

32 

4.748 

4.730 

4.730 

4.729 

4.723 

33 

4.739 

4.732 

4.722 

4.721 

4.716 

3^ 

4.734 

4.731 

4.719 

4.718 

4.711 

35 

4.729 

4.733 

4.720 

4.720 

4.713 

3  6 

4.727 

4.732 

4.716 

4.716 

4.707 

37 

4.729 

4.717 

4.722 

4.722 

4.715 

38 

4.725 

4.710 

4.710 

4.710 

4.701 

39 

4.717 

4.715 

4.727 

4.728 

4.723 

4.769 

4.667 

4.663 

4.664 

4.661 

41 

4.698 

4.472 

4.738 

4.738 

4.728 

42 

3.948 

3.908 

3.788 

3.788 

3.770 

43 

2.973 

3.087 

3.030 

3.030 

3.025 

44 

2.863 

2.887 

2.884 

2.884 

2.883 

45 

2.857 

2.861 

2.861 

2.861 

2.861 
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Table  4  Velocity 


Cell 

e=a 

x=o 

e=i 

1=0.5 

1 

0.780 

0.775 

2 

0.793 

0.788 

3 

0.819 

0.813 

4 

O.858 

0.851 

5 

0.909 

0.900 

6 

0.972 

0.961 

7 

1.044 

1.032 

8 

1.125 

1.106 

9 

1.210 

1.195 

10 

I.298 

1.282 

U 

1.382 

1.365 

12 

1.456 

1.440 

13 

1.511 

1.498 

14 

1.542 

1.534 

15 

1.551 

1.547 

16 

1.546 

1.546 

17 

1.541 

1-543 

18 

1.540 

1.540 

19 

1-539 

1.538 

20 

1.539 

1.536 

21 

1.538 

1.535 

22 

1-538 

1.534 

23 

1-538 

1-533 

24 

1.538 

1.533 

25 

1.538 

1.532 

26 

1.538 

1.532 

27 

1.538 

1.532 

28 

1.538 

1-532 
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Table  4  (Cont. ) 


29 

1.538 

1.532 

30 

1.538 

1.532 

31 

1.538 

1.532 

32 

1-538 

1.532 

33 

1-538 

1.532 

34 

1.538 

1.532 

35 

1.538 

1.532 

36 

1.538 

1.531 

37 

1.538 

1.529 

38 

1.538 

1.521 

39 

1-537 

1.489 

4o 

1.552 

1.382 

4i 

1.177 

1.079 

4 2 

0.262 

0.581 

43 

0.029 

0.213 

44 

0.003 

0.064 

45 

0.000 

0.018 

46 

0.000 

0.005 

47 

0.000 

0.001 

48 

0.000 

0.000 

49 

0.000 

0.000 

50 

0.000 

0.000 
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EXACT  SOLUTION 


Figo  1.  Velocity  profile  at  time  t  =  13.48  for  problem  1, 
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Fig.  2.  Density  profile  at  time  t  =  13.48  for  problem  1. 
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Fig*  3*  Velocity  profile  of  the  shock  front  at  t  =  13*46  for  various 
choices  of  X. 


